Nonlinear  Fluid  Flow/Surfactant/Interface  Dynamics 


Antony  N.  Boris 

Department  of  Chemical  Engineering 
University  of  Delaware 
Newark,  DE  19716 

phone:  (302)  831-8018  fax:  (302)  831-1048  e-mail  Address:  beris@che.udel.edu 

Award#:  N00014-98- 1-0008 


LONG-TERM  GOALS 

The  goal  of  this  research  is  the  development  of  a  theoretical  quantitative  understanding  of  the  dynamics 
of  high  Reynolds  number  free-surface  flows  under  large  free-surface  deformations  in  the  presence  of 
one  or  multiple  surfactants.  The  long-term  goal  of  this  research  is  to  understand  the  effect  of  the 
surfactants  on  (1)  turbulence,  (2)  waves,  (3)  slick  formation,  and  (4)  mass  transfer  through  the  free 
surface. 

OBJECTIVES 

The  objectives  of  the  proposed  research  are  several.  Eirst,  we  are  aiming  to  develop  an  algorithm  for 
direct  numerical  simulation  of  free-surface  flows  with  surfactants.  The  scheme  will  be  capable  of 
simulating  two-dimensional  flows  with  high  Reynolds  numbers,  accommodating  nonlinear  free-surface 
deformation  and  the  presence  of  one  or  two  interacting  surfactants.  Second,  we  want  to  use  this  tool  to 
evaluate  alternative  approaches  in  modeling  surfactant  behavior.  In  particular,  through  comparison 
with  experimental  data,  we  want  to  investigate  the  extent  of  microscopic  detail  needed  to  be  present  in 
the  surfactant  model  in  order  to  enable  an  at  least  qualitative  description  of  the  dynamic  surfactant 
behavior.  We  want  to  focus  on  improving  the  understanding  of  the  flow/free-surface/surfactant 
interactions  through  a  direct  simulation  of  the  nonlinear  interactions  of  a  surfactant-contaminated  free 
surface  with  a  vortex  with  or  without  the  presence  of  waves. 

APPROACH 

The  approach  followed  for  the  accurate  and  computationally  efficient  evaluation  of  the  flow,  free- 
surface  deformation  and  surfactant  concentration  involves  to  development  of  a  direct  numerical 
simulation  capability  based  on  a  fully  spectral  spatial  discretization  and  a  conformal  mapping  of  the 
flow  domain  into  a  regular  parallelepiped.  The  conformal  mapping  is  key  in  order  to  preserve  the 
computational  efficiency  of  the  proposed  numerical  technique  since  it  allows  the  use  of 
computationally  efficient  fast  Poisson  algorithms  for  the  solution  of  the  Poisson  and  Helmholtz 
problems  generated  from  the  time-integration  of  the  flow  and  surfactant  concentration  evolution 
equations.  The  numerical  time-integration  approach  can  be  either  a  mixed  implicit/explicit  scheme,  or 
a  fully  implicit  scheme.  These  approaches  have  led  to  very  efficient  numerical  schemes  for  stationary 
boundary  problems.  The  mixed  integration  scheme  has  recently  been  extended  to  three  dimensions. 

Eor  free-surface  problems,  we  have  developed  a  fully  implicit  scheme  for  simulating  fully  nonlinear 
free-surface  problems.  We  are  currently  incorporating  the  capability  to  handle  surfactants,  which  is 
now  straightforward.  In  addition  to  this  work,  we  have  gained  extensive  experience  in  parallel 
computing  and  have  the  capability  of  performing  our  simulations  in  powerful  parallel  computing 
environments  through  very  efficient  implementations  of  our  algorithms,  which  are  linearly  scalable 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 
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with  the  number  of  processors.  Thus,  we  are  in  a  unique  position  to  perform  efficient  simulations  of 
free-surface  flows  with  surfactants,  accommodating  fully  nonlinear  free-surface  deformations,  with 
very  high  accuracy. 

WORK  COMPLETED 


During  the  past  year,  we  extended  significantly  our  algorithm  for  flow  in  channels  with  deformed 
boundaries  by  improving  its  range  of  applicability  in  three  directions,  retaining  its  almost  linear 
scalability  (the  number  of  operations  required  are  order  Nlog2N,  where  N  is  the  number  of  unknowns): 

•  We  successfully  extended  the  solid  boundary  algorithm  to  incorporate  a  third  dimension. 

•  We  implemented  two  different  schemes  for  fully  implicit  time  integration. 

•  We  successfully  extended  our  two-dimensional  algorithm  to  accommodate  a  free  surface. 


Specifically,  for  our  three-dimensional  algorithm  we  utilized  the  mixed  explicit/implicit  integration 
scheme  discussed  in  a  recent  publication  in  the  Journal  of  Computational  Physics  dealing  with  the 
two-dimensional  case  (Dimitropoulos  et  ah,  1998a),  extending  it  through  the  use  of  a  modified 
preconditioner  used  in  the  iterative  solver.  The  main  difficulty  in  the  three  dimensional  case  is  the 
efficient  implementation  of  the  influence  matrix  method  for  satisfying  the  divergence-free  condition. 
For  the  type  of  mapping  used  in  our  case  (pseudoconformal  mapping),  the  divergence  now  has  a 
function  (a  scale  factor)  multiplying  the  third  term  corresponding  to  the  spanwise  direction  (z),  which 
leads  to  the  of  coupling  the  Fourier  modes  and  to  a  non-separable  equation: 
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where  (^,  p,  z)  are  the  pseudoconformal  coordinates  and  M  is  the  constant  ratio  of  the  scale  factors 
along  ^  and  p.  The  scale  factor  along  z  is  equal  to  unity.  To  get  around  this  obstacle,  we  decided  to 
satisfy  the  divergence  iteratively  within  the  preconditioner,  utilizing  a  similar  Concus  and  Golub  type 
transformation  as  for  the  Navier-Stokes  equations.  This  in  principle  leads  to  an  convergent  algorithm 
for  flows  which  do  not  have  dominating  flow-field  variations  in  the  spanwise  direction  as  one  can  see 
from  our  analysis  of  the  basic  scheme  in  our  previous  work  (Dimitropoulos  and  Beris,  1997; 
Dimitropoulos  et  ah,  1998a).  This  algorithm  is  currently  being  validated  for  flows  with  three- 
dimensional  flow-fields  in  a  parallel  supercomputer.  Additionally,  we  now  perform  our  simulations  on 
parallel  machines  using  MPI,  where,  as  expected,  we  observe  almost  linear  scalability  of  performance 
with  the  number  of  processors,  consistent  with  the  initial  design  objectives  for  our  algorithms. 


For  the  free-surface  problem,  we  have  observed  that  stability  necessitates  the  use  of  fully  implicit 
algorithms  when  dealing  with  nonlinear  boundary  conditions.  We  have  implemented  implicit  time- 
integration  in  two  different  ways.  First,  we  used  second  order  Adams-Moulton  time-integration  for  the 
nonlinear  terms  in  the  Navier-Stokes  equation,  which  leads  to  nonlinear  equations  to  be  solved  at  each 
time  step.  We  proceed  to  solve  this  system  of  equations  by  combining  our  algorithm  with  an  external 
Newton-Raphson  scheme.  This  is  especially  suited  for  use  in  problems  with  large  contributions  from 
the  nonlinear  terms.  Our  procedure  consists  of  linearizing  our  equations  and  attaining  the  form: 

V"(dv)-  ^^(dv)-  2  Re  y{dp)^  RHS  (t  <  ?"*')+  RHS  v,  p)+  RHS  (t"*^-S\,Sp} 

where  the  right-hand  sides  (RHS)  can  be  thought  of  as  a  contribution  belonging  to  the  previous  steps,  a 
contribution  from  the  current  time-step  which  is  updated  at  every  Newton’s  iteration  and  a  contribution 
which  is  updated  continuously  within  each  Newton’s  iteration  as  we  solve  for  the  corrections  5q,  where 


q  is  any  variable  through  an  iterative  algorithm,  the  preconditioned  conjugate  gradient  algorithm 
described  in  our  earlier  work.  This  method  has  no  problems  in  satisfying  the  divergence-free  condition 
since  each  iterate  is  found  to  be  divergence-free  by  the  same  influence  matrix  algorithm  as  before.  The 
computational  cost  per  time- step  of  this  scheme  compared  to  our  previous  algorithm  for  two- 
dimensional  problems  is  proportional  to  the  number  of  Newton-Raphson  iterations  used,  typically 
around  3-5  when  a  decent  and  divergence-free  initial  guess  is  supplied. 

One  can  also  implement  an  implicit  approach  using  a  predictor-corrector  scheme,  using  an  Adams- 
Bashforth  scheme  for  the  non-linear  terms  in  the  predictor  step  and  an  Adams -Moulton  scheme  for  the 
corrector  step.  This  algorithm  is  only  2  times  more  expensive  than  the  mixed  scheme  and  has  the 
advantage  that  the  initial  guess  need  not  satisfy  the  continuity  equation  very  well,  allowing  for  more 
freedom.  Of  course,  the  main  disadvantage  compared  to  the  full  Adams-Moulton  scheme  using  the 
external  Newton-Raphson  scheme  is  that  it  is  not  as  stable.  However,  it  is  adequately  stable  for  use 
even  in  free-surface  problems.  The  concept  behind  the  algorithm  is  simple.  One  recalculates  the 
nonlinear  terms  at  the  current  time  step  after  the  predictor  step  is  finished  and  resolves  with  a  modified 
right-hand  side  corresponding  to  these  terms.  In  fact,  the  corrector  step  converges  faster  than  the 
predictor  step,  since  the  iterative  algorithm  starts  with  a  better  initial  guess.  Therefore,  the  algorithm  is 
actually  less  than  a  factor  of  2  slower  than  that  using  a  mixed  time-integration  scheme. 

The  solution  of  the  free-surface  problem  was  initially  attempted  through  the  use  of  the  mixed 
algorithm.  However,  stability  considerations  made  us  abandon  that  scheme  in  favor  of  a  predictor- 
corrector  scheme,  which  can  handle  the  problem  in  the  least  intrusive  manner  compared  to  the  solid 
boundary  code.  Starting  from  our  code  for  moving  solid  boundaries  (see  Dimitropoulos  et  al,  1998a), 
we  can  implement  the  kinematics  condition  along  with  the  orthogonality  constraint  in  the  algorithm 
that  computes  the  time-derivatives  of  the  mapping.  This  is  achieved  by  solving  iteratively  two 
separable  Poisson  equations,  using  a  similar  method  to  that  in  Dimitropoulos  et  al.  (1998a),  with  the 
following  (coupled)  boundary  conditions: 
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Once  the  time-derivatives  of  the  mapping  are  known,  we  can  then  find  the  actual  mesh  by  integrating 
in  time.  To  enforce  orthogonality  over  the  time-integration  error,  we  fit  the  top  boundary  calculated 
through  the  time-integration  to  a  Fourier  series  and  then  use  this  Fourier  series  as  afunctional 
representation  in  the  mapping  algorithm  for  a  given  surface  (described  in  Dimitropoulos  et  al.,  1998a). 
Vke  then  proceed  to  solve  for  the  flow  problem  with  a  known  mesh  at  the  new  time-step.  This  is  made 
robust  through  a  predictor-corrector  scheme.  Initially,  we  use  explicit  integration  for  the  mesh.  We 
solve  for  the  velocities  and  then  recalculate  the  time-derivatives  of  the  mesh  with  this  new  velocity  and 
use  implicit  integration  for  the  mesh.  We  proceed  to  solve  again  for  the  velocity  field.  Once  this  first 
corrector  step  is  completed,  we  repeat  it  until  the  values  we  obtain  do  not  change.  This  is  necessary 
since  we  have  decoupled  the  solution  of  the  mesh  from  the  solution  of  the  flow  problem.  This  would 
not  be  possible  using  the  fully  implicit  scheme  with  the  Newton-Raphson  iteration.  This  approach 
allows  us  to  implement  the  full  nonlinear  free-surface  boundary  conditions  without  needing  to  change 
substantially  the  preconditioned  conjugate  gradient  solver.  The  free-surface  boundary  conditions  as 
expressed  in  the  pseudoconformal  coordinate  system  utilized  in  this  work  are: 
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tangential  stress, 


where  Ca  is  the  capillary  number  and  H  is  the  mean  curvature  of  the  interface.  We  implement  these 
boundary  conditions  in  the  preconditioner,  along  with  the  divergence-free  condition  after  transforming 
them  into  the  following  set  of  equations: 
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The  constant  K  is  simply  the  average  of  the  square  of  the  scale  factor,  RBCl,  RBC2,  RBC3  are  the 
residuals  of  the  boundary  conditions  calculated  in  the  conjugate  gradient  level  and  i  is  the  index  of  the 
Concus  and  Golub  type  iteration  performed  in  the  preconditioner.  This  scheme  differs  from  that  for  the 
solid  boundary  case  in  that  the  boundary  conditions  change  at  every  preconditioner  iteration. 


RESULTS 


We  present  briefly  a  validation  result  for  the  free-surface  simulation  algorithm.  We  simulated  a  wave 
of  amplitude  equal  to  10'^  in  a  channel  2  units  deep  and  1  unit  wide.  The  Reynolds  number  used  was 
10**  (based  on  half-depth  and  gravity),  whereas  the  corresponding  capillary  number  for  water  was  10  ‘. 
The  time-step  was  equal  to  lO  l  This  corresponds  to  a  period  of  2.55  time  units.  Our  results  agree  well 
with  the  theoretical  results  of  Levich  (1960)  for  damping  of  waves  with  small  viscosity.  Figure  1 
shows  the  time-evolution  of  the  interface  shape.  Detailed  results  for  all  the  work  discussed  in  this 
report  will  be  available  in  a  series  of  future  publications. 

IMPACT/APPLICATIONS 


We  have  extended  our  computationally  efficient  numerical  scheme  for  solving  spectrally  two- 
dimensional,  time-dependent  flow  problems  in  moderately  complex  geometries  to  accommodate 
implicit  time-integration,  a  third  dimension  for  the  solid  boundary  case  and  a  free-surface  for  the  two- 
dimensional  case.  Our  methods  exhibit  almost  linear  scalability  and  exponential  convergence  with 
mesh  refinement.  They  use  an  orthogonal  mapping  algorithm,  an  efficient  and  robust  iterative  solver 
incorporating  the  influence  matrix  method  for  satisfying  the  incompressibility  condition.  In  addition 
the  algorithms  are  parallelizable  in  a  straightforward  fashion  and  exhibit  linear  scalability  of  their 
performance  with  the  number  of  processors  used.  They  provide  unique  new  computational  tools  for 
the  calculation  of  complex  multidimensional  and  time-dependent  flows,  such  as  free-surface  flows  with 
surfactants,  turbulent  flow  over  irregular  (albeit  smooth)  boundaries,  blood  flow  in  arteries  and 
polymer  resin  flow  in  model  porous  media,  which  is  of  importance  to  composite  manufacturing.  We 
are  currently  extending  the  free-surface  code  for  simulation  of  free-surface  flows  with  surfactants, 
where  we  solve  the  fully  coupled  nonlinear  problem  of  fluid  flow  and  mass  transport  of  the  surfactant. 
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Figure  1:  Evolution  in  time  of  the  interface  shape  of  a  small-amplitude, 
large-wavelength  water  wave. 


TRANSITIONS 

The  main  development  of  this  work  is  a  quite  general  in  applieability  algorithm,  whieh  extends  the  use 
of  full  spectral  methods  to  flows  within  moderately  complex  geometries.  The  tools  developed  in  this 
work  i.e.  the  pseudoconformal  mapping  and  the  efficient  spectral  preconditioner  are  also  applicable  to 
other  computational  fluid  mechanics  and  transport  phenomena  problems.  All  of  these  are  expected  to 
find  significant  use  within  the  computational  fluid  dynamics  community,  since  their  possible 
applications  span  a  quite  wide  area. 

RELATED  PROJECTS 

The  core  of  the  spectral  numerical  method,  after  suitable  modification,  has  been  successfully  used  in 
direct  numerical  simulations  (DNS)  of  multidimensional  and  time-dependent  viscoelastic  flows  in  a 
joint  University  of  Delaware  -  NRL  effort  for  the  investigation  of  polymer-induced  turbulent  drag 
reduction.  We  have  gained  significant  experience  in  parallel  computing  from  this  project,  since  we 
have  performed  extensive  simulations  on  the  CRAY  T3D,  T3E  systems  at  the  Pittsburgh 
Supercomputing  Center  (PSC)  and  the  CRAY-Origin2000  at  National  Center  for  Supercomputing 
Applications  (NCSA).  This  work  has  led  to  a  publication  discussing  our  first  calculations  in  Physics  of 
Fluids  (Sureshkumar  et  ah,  1997)  and  another  in  the  Journal  of  Non-Newtonian  Fluid  Mechanics 
(Dimitropoulos  et  al,  1998b),  where  the  effect  of  rheology  in  polymer- induced  drag  reduction  is 
examined.  The  transfer  of  the  codes  developed  in  the  current  project  from  the  single-processor  to 
multiple-processor  environments  has  been  greatly  facilitated  from  the  above-mentioned  experience. 
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